Graphene on Ir(lll) surface: Prom van der Waals to strong bonding 
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We calculate the properties of a graphene monolayer on the Ir(lll) surface, using the model 
in which the periodicities of the two structures are assumed equal, instead of the observed slight 
mismatch which leads to a large superperiodic unit cell. We use the Density Functional Theory 
approach supplemented by the recently developed vdW-DF nonlocal correlation functional. The 
latter is essential for treating the van der Waals interaction, which is crucial for the adsorption 
distances and energies of the rather weakly bound graphene. When additional iridium atoms are 
put on top of graphene, the electronic structure of C atoms acquires the sp 3 character and strong 
bonds with the iridium atoms are formed. We discuss the validity of the approximations used, and 
the relevance for other graphene-metal systems. 
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I. INTRODUCTION 

Graphene is a one-atom thick two-dimensional struc- 
ture of carbon atoms arranged in a honeycomb lat- 
tice. It is (conceptually at least) at the origin of all 
other graphitic forms ji including the three-dimensional 
graphite, one-dimensional carbon nanotubes and zero- 
dimensional fullerenes. The planar geometry and the ex- 
ceptional strength of graphene^ are due to the sp 2 bonds 
between atoms. Single-layer graphene has been obtained 
by micromechanical cleavage of graphite and by growth 
on SiC and metal surfaces. The recent large increase of 
interest in graphene is due both to the theoretical im- 
plications of its unique electronic properties and to its 
potential applicability, in particular as a novel material 
for electronics. 

As the building block of graphite and as the adsorbate 
on many surfaces, graphene bonds to its surroundings 
only weakly, and the character of the bonding is largely 
van der Waals (vdW) . Occasionally, stronger bonding oc- 
curs without destroying the geometry of the graphene lat- 
tice^ for example on Ni(lll)£ and Ru(0001)£~— surfaces. 
Graphene on Ir(lll) is an example where, depending 
on conditions, both kinds of bonding can occur. Mono- 
layer graphene is vdW physisorbed, and the characteris- 
tic graphene lattice can be clearly seen in STM images j& 
while with additional Ir clusters on top the carbon-metal 
bonds become stronger. 

Large-cell Density Functional Theory (DFT) calcula- 
tions of graphene on Ir(lll) using PBE GGA (Ref. i) 
and LDA (Ref. 0) have been performed. However, exper- 
iments and calculations reveal a subtle interplay between 
van der Waals bonding and stronger electronic interac- 
tion with the substrate. The vdW interaction is not prop- 
erly described in the standard local (LDA) and semilocal 
(GGA) DFT functionals, which provides a serious obsta- 
cle to the complete understanding of the nature of the 
bonding. 

In this paper we apply the recently developed extension 
to the DFT, which replaces the semi-local (i.e. depend- 



ing upon the gradient of the electronic density) correla- 
tion term with a fully non-local one (depending upon the 
electronic densities at different points in space), which 
can describe the van der Waals forces even between two 
fragments of matter with non-overlapping electronic den- 
sities. Due to computational complexity we had to aban- 
don the large supercell which aims to describe more re- 
alistically the graphene and Ir(lll) surface with their 
slightly different atomic periodicities, and opt for an ap- 
proximate description by a smaller commensurate unit 
cell. While the quantitative accuracy of the results suf- 
fers (but we argue that it is a quite limited and controlled 
problem), the approximation makes the transition from 
weak to strong bonding easier to analyze and understand. 



II. GRAPHENE BINDING IN GRAPHITE 




FIG. 1: The structure of graphite, side and top view. The 
layers are stacked in AB order, so that half of the C atoms lie 
in chains along the direction perpendicular to the graphene 
planes, while the other half alternates in the other two high- 
symmetry positions. The periodicity in the perpendicular di- 
rection is c, i.e. the interlayer distance is c/2. For clarity, c 
has been exaggerated by about a factor of 3. 
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FIG. 2: Binding energies of graphite AB structure as a 
function of interlayer distance: Comparison of various DFT 
codes, plane- wave based dacapo— and abinit,— and real- 
space gpaw.— Local LDA and semilocal PBE GGA results 
are shown. The curve labeled vdW-DF is the energy when 
the PBE GGA correlation has been replaced by a fully non- 
local correlation 1 - using the JuNoLo numerical code^. 

We have first applied our methods to graphite, i.e. 
graphene sheets arranged in the most stable AB stack- 
ing, shown in Fig. [TJ This is a much studied system with 
good experimental data and calculated values of struc- 
tural and energetic parameters. It has the same complex- 
ities of having both the weak vdW and strong chemical 
bonds as our principal subject of interest, graphene on 
Ir(lll). A single sheet of graphene presents no difficul- 
ties for the standard DFT GGA approach, giving the C-C 
distance of 1.42 A, corresponding to a lattice constant of 
2.46 A, in agreement with the experiment. Next, we per- 
formed the Density Functional calculations of stacks of 
graphene sheets using several flavors of LDA and GGA, 
implemented in various numerical codes. The calculated 
cohesive energies are shown in Fig. [21 as a function of 
the interlayer separation. After that we used the nonlo- 
cal vdW-DF functional for the correlation! 13 ' 15 

The pure DFT results agree well for all programs used, 
and are even quite insensitive upon the (lack of) full 
self-consistency. For example, dacapo calculations use 
the PW91 GGA functional, but the dacapo LDA curve 
shown in the figure, obtained evaluating the LDA func- 
tional on the electron density calculated with PW91, 
agrees quite well with the fully selfconsistcnt abinit LDA 
results. GGA calculations give little or no bonding, 
and LDA gives (apparently) reasonable bonding ener- 
gies and distances, comparable to experimental values. 
The reason for the failure of GGA is intuitively clear: 
The semilocal gradient approximation cannot describe 
well the inherently non-local van der Waals interaction, 
which exists even between subsystems with completely 
non-overlapping electronic densities. The apparent suc- 
cess of LDA is somewhat perplexing, since it is even more 



local than the more advanced GGA, the latter indeed be- 
ing more successful when it comes to chemically bound 
systems. There are strong indications that the agreement 
of the LDA results is largely fortuitous, as discussed fur- 
ther on. 

We have further investigated the problem by apply- 
ing the vdW-DF theory^ which is at present the most 
promising approach for treating the non-local correlation. 
It consists in replacing the semilocal (gradient) part of 
the GGA correlation functional by a fully non-local term, 
which still depends only upon the electronic density, in 
the true spirit of the Density Functional theory. We have 
applied vdW-DF as implemented in the JuNoLo codei^ 
in a post-processing approach, i.e. we used the electron 
densities obtained in the standard GGA calculation to 
evaluate the non-local vdW-DF correlation, and inserted 
it into the total energy instead of the semilocal correla- 
tion. This approach is, of course, not fully selfconsistcnt, 
since the DFT potential and the Kohn-Sham wavefunc- 
tions, and hence the electron density can depend upon 
the details of the correlation functional used. However, 
a recent selfconsistcnt implementation of the vdW-DF 
correlation functional shows that the differences are neg- 
ligible^ We have therefore relied on the post-processing 
approach which is less time consuming and avoids any in- 
tervention in the code of the DFT programs. Changing 
the correlation contribution changes the total energy, and 
therefore the forces acting on the atoms as well. How- 
ever, all atomic configurations which we consider here 
have high symmetry, where we can sweep the interest- 
ing range of interlayer separation "by hand" in order to 
find the optimum configuration. The lack of fully selfcon- 
sistent atomic relaxations inherent to such an approach 
is not a major problem, as discussed later on. We fur- 
thermore note that we have not followed the suggestion 
put forward by the authors of the vdW-DF theory to 
use the revPBE exchange functional^ 3 - and have instead 
continued using the PBE exchange. Although revPBE 
exchange seems to compensate for too large binding en- 
ergies for several van der Waals bound systems, it gives 
worse equilibrium distances, and the same improvement 
does not seem to occur in cases of strong bonding. 

The vdW-DF results shown by a thick line in Fig. [2] 
are qualitatively similar to the LDA results, but with 
some important differences. The equilibrium distance at 
around 3.6 A is larger than the LDA result, as is the 
binding energy of 30.8 meV per carbon atom. The vdW- 
DF attractive potential has clearly a longer range than 
LDA, which reflects the long-range nature of the van 
der Waals attraction and reveals the fortuitous character 
of the agreement with LDA around the minimum. The 
vdW-DF values for the interlayer separation and inter- 
action energy agree well with recent experimental results 
and theoretical calculations (See Ref. [l7| and references 
therein). Recently, it has been found that the behavior 
of the nonretarded van der Waals interaction between 
nonovcrlapping anisotropic nanostructures that have a 
zero electronic energy gap should be different than pre- 
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dieted from the the usual sum of R~ 6 contributions^ 
but this is probably relevant only in the extreme asymp- 
totic regime. 

We conclude that the inclusion of the van der Waals 
interaction is essential to reproduce physical properties 
of graphite, and that the vdW-DF approach successfully 
treats all aspects of the graphene binding in graphite. 



III. STRUCTURE OF GRAPHENE ON Ir(lll) 
SURFACE 




FIG. 3: Moire superstructure of 10 x 10 graphene on 9 x 9 
Ir(lll) unit cell. The C atoms are (approximately) above the 
first and third layer Ir atoms within the circle labeled 1-3, 
above the first and second layer in 1-2, and above the second 
and third layer in 2-3.— 

Graphene monolayers of high structural quality, ex- 
tending over tens of nanometers and even up to microm- 
eter size, orientationally well aligned with the substrate, 
have recently been obtained by hydrocarbon decomposi- 
tion on Ir(lll)4 The lattice constants of graphene and 
the Ir(lll) surface differ at room temperature by around 
10%, and the STM micrographs clearly show the moire 
pattern due to the lattice mismatch. In Fig. [3] the su- 
percell with a 10 x 10 graphene lattice on top of 9 x 9 
structure of iridium atoms is shown. A further intrigu- 
ing feature is observed when additional iridium atoms are 
adsorbed on top of the graphene. STM images show that 
the adatoms form regular arrays on clusters, selectively 
bound to certain regions of the moire pattern4^ 

Density functional calculations employing the PBE 
GGA functional^ and LDA functional^, on a supercell 
similar to the one in Fig. [3] have been performed. It has 
been found thet the PBE GGA functional gives almost 
no bonding of graphene monolayer on Ir(lll) (Refs. HI 
Erratum. l20t). while the LDA functional gives reasonable 
results for bonding energies and interatomic distances, 
both without and with additional clusters on top^ These 
results are qualitatively reminiscent of our results for 
graphite, i.e. we again see an apparent success of the 
quite basic LDA, which signals that a more detailed in- 
vestigation is necessary. 

Our strategy is similar to the approach used for 
graphite in Section [TTJ We start with the standard DFT 



and later investigate the effects of the non-local correla- 
tion. We do not use the realistic large unit cell shown 
in Fig. 02 but instead we compress the Ir(lll) substrate 
in the surface plane (coordinates x, y) so that it matches 
the lattice constant of graphene. By changing the phase 
of the carbon atoms with respect to the underlying lat- 
tice of Ir atoms, we are able to simulate (approximately) 
any point in the supercell in Fig. [3J We have done most 
calculations for the region labeled 1-3, which both ex- 
periment and calculations show to be the most strongly 
bonding, with only a few checks of the other regions. 

Calculations of commensurate graphene-metal surface 
systems has been done for several metal surfaces, ei- 
ther by adjusting the substrate lattice constant^! or the 
graphene lattice constant 2 ^. We shall discuss these cal- 
culations into more detail later on. 

The mismatch of the lattice constant of graphene 
(2.46 A) and that of the Ir(lll) surface (2.73 A, corre- 
sponding to the conventional fee lattice constant clq = 
3.86 A) is around 10%, clearly larger than those in 
Ref. [2l|, which arc in the range of 0.8-3.8%. In order 
to minimize possible artefacts due to the squeezing of 
the iridium substrate to fit the graphene lattice we opti- 
mized the lattice constant of the iridium substrate in the 
z direction. To that end, we performed calculations of 
iridium bulk with compressed (111) planes and allowed 
it to relax in the perpendicular direction. The lattice 
constant in the z direction increased by about 10% to 
4.24 A, and we used this lower- symmetry iridium lattice, 
compressed in x — y directions and expanded in z, as the 
substrate in our calculations. From the point of view of 
quantitative accuracy, a calculation using a large super- 
cell and "natural" iridium substrate would be preferred, 
but our approach enables a clear insight into the bonding 
properties, which would be at risk to remain buried and 
hard to see in a more realistic large calculation. 

Another important aspect of graphene interaction with 
the Ir(lll) surface can be inferred from the band struc- 
ture of Ir(lll) along high-symmetry directions of the 
surface Brillouin zone. Our calculations based on the 
DFT Kohn-Sham eigenstates^ as well as ARPES exper- 
iment o 23 ' 24 show that there is an energy gap around the 
K point of the surface Brillouin zone, extending from 
just below the Fermi level down to almost 1.5 eV bind- 
ing energy. The vertex of the "Dirac cone" of the ir bands 
of graphene adsorbed on Ir(lll) lies entirely within this 
gap^l The weak interaction of graphene with the iridium 
substrate can be attributed to this mismatch of the elec- 
tronic states, since the unsaturated ir bands of the Dirac 
cones don't have any substrate states with the same mo- 
mentum k and energy E to hybridize with. We have 
also checked the band structure of our compressed irid- 
ium surface and found that the band gap around the K 
point is still present and has a similar shape, which im- 
plies that the weak character of the graphene interaction 
with Ir(lll) will not be much affected by the use of the 
compressed substrate. 

In our DFT calculations we use a three-layer Ir(lll) 
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slab with the adjusted lattice constants in the x — y and 
the z directions as explained above. It would have been 
easy to use a thicker substrate, but this would add no 
further quantitative accuracy, considering other simplifi- 
cations and approximations used. 



IV. DFT CALCULATIONS OF GRAPHENE ON 
Ir(lll) AND Ir-GRAPHENE-Ir SANDWICHES 

For our calculations we have chosen four characteristic 
structures of graphene on iridium, shown in Fig. [4j The 
structures are periodic in the x — y plane and extend to 
infinity, but for clarity we show only a small symmetric 
cluster of atoms for each one. There are two atoms in the 
unit cell of graphene, which in the following we denote 
by Ca and Cb, as illustrated in Fig. [4] (d). We denote 
the iridium atoms in the first substrate layer by Irs, and 
the atoms in the first overlayer as Iro^ 

The structures were chosen so that they illustrate Ir-C 
bonds of various character, with overall bonding strength 
increasing from structure (a) to structure (d). The irid- 
ium substrate is modeled by three atomic layers in all 
cases. The structures are: (a) Monolayer graphene on 
Ir(lll) with Ca above Irs and Cb above third layer Ir, 
which illustrates the most stable regions of the moire pat- 
tern of a monolayer graphene on Ir(lll); (b) Graphene 
monolayer as in (a) , but with a single overlayer of iridium 
atoms, with Iro located above the centers of the hexag- 
onal rings of the graphene; (c) Three additional layers 
of iridium, with Iro above Cb; (d) A single Ir overlayer, 
with atoms in the same positions as in the first layer in 
(c). In (c) and (d) there is one iridium atom below Ca 
and another one above Cb , which models the geometries 
of the stable iridium clusters on top of graphene observed 
in the experiment. 
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(a) 

FIG. 4: The four structures considered in the paper. For 
clarity, only a few atoms from each atomic layer are shown in 
structures (b)-(d). 



A. Standard DFT only 

We first calculated the dependence of the interaction 
energy upon Ir-graphene separation for structures (a)- 
(d), using the standard LDA and GGA PBE functionals. 




FIG. 5: DFT (LDA and PBE GGA) energies for graphene 
on Ir, structure (a) in Fig. [4l and Ir-graphene- Ir sandwiches, 
structures (b)-(d) in Fig. 2] In this and the following graphs 
the energies are given per unit cell, i.e. two C atoms and one 
Ir atom in each iridium layer. 



The results are shown in Fig. [5] Two comments are in or- 
der. For structure (a), the energy scale in Fig.[5]is smaller 
by a factor of two, since the graphene has iridium atoms 
only on one side, and the interaction (in particularly the 
repulsion at small distances) is expected to scale with 
the number of neighboring atomic planes. Secondly, for 
the sandwich structures (b)-(d) the graphene layer was 
at the beginning of the calculations placed symmetrically 
between the nearest Ir layers, each of them at a distance 
of zi r _g, but was allowed to relax in the course of the cal- 
culation. At small separations zi r - g the graphene buck- 
les, with Ca atom moving towards Irg and Cb towards 
the overlayer. This is not a major problem, since the 
two atoms move by almost the same amount in oppo- 
site direction, even in the case of the less symmetrical 
structure (b), so that zi r -g still measures the ^-averaged 
position of the graphene plane. There is no buckling 
at zir-g separations larger than say 3 A. Since the re- 
laxation was done according to the forces calculated in 
GGA functional, which does not develop any apprecia- 
ble attractive potential well at these distances, there was 
no danger that the graphene layer would be attracted to 
the iridium atoms on one side. Iridium atoms were not 
relaxed. 

The results for all structures almost coincide for zi r - g 
larger than around 3.3 A, for both GGA and LDA, i.e. 
the interaction at physisorption distances does not show 
a large corrugation along the x — y coordinates. The 
sandwich structures (b)-(d) show the tendency to form 
a strong bond at small graphenc-Ir distances. The min- 
imum is around 2.3 A, and is quite sensitive on the rel- 
ative position of the atoms in the x — y plane. Thus in 
the unfavorable structure (b), where the Ir atoms in the 
additional layer do not lie directly above the C atoms, 
there is only a kink in the interaction energy at about 
2.3 A, hinting that there is a tendency towards strong 
chemical bonding. The structures (c) and (d) develop a 
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distinct potential well around that distance, which is not 
deep enough in the PBE GGA functional calculation, but 
with the LDA functional it becomes the stable configu- 
ration with more than 0.5 eV binding energy. Note that 
the quantity zi r - g measures the distance to the average 
z coordinate of the graphene layer, and since there is a 
buckling of about 0.2 A of C atoms towards the nearest 
Ir atom, the Ir-C distance is actually around 2.1 A, as 
discussed more into details later on. 

These results prompt us to reevaluate even the stan- 
dard DFT calculations in the region of strong bonding 
at small Ir-C distances, where the graphene layer sig- 
nificantly changes its electronic character. Up to now 
we have consistently used the lattice constant of free 
graphene, assuming that it is optimal for the bound sys- 
tem too (and we went to the trouble of compressing the 
Ir layers accordingly) . This may not be true in the region 
of strong bonding, and we first check this. 

To that end, we performed standard DFT calculations 
of structures (c) and (d) with the lattice constant in 
the x — y plane slightly expanded from the value of free 
graphene (3.478 A, C-C distance 1.42 A) and accordingly 
reduced in the z direction, and checked whether there 
was any improvement of the total energy. In order to 
evaluate the bonding energy we also had to calculate the 
energy of separated Ir and graphene slabs (corresponding 
to Zir- g — > oo) for each value of the expanded lattice con- 
stant, and subtract it from the energy of the interacting 
system. In Fig. [7] we show the results for structures (c) 
and (d). The unconnected circles are the non-optimized 
results for the two structures taken from Fig. [5j while 
the connected circles are the best results obtained by ex- 
panding the lattice. We see that the energy improves 
significantly in the region of strong bonding, where the 
optimum lattice constant at the position of the minimum, 
Zir- g ~ 2.3 A, increases from the free graphene value of 
3.478 A to 3.65 A for structure (c) and to 3.72 A for 
structure (d). The fact that in the region of weak bind- 
ing, for Zir-g > 3 A, there is no improvement of energy 
and the optimal lattice constant remains at the value 
of free graphene (i.e. the large graphene stiffness dom- 
inates the energy balance) indicates that the procedure 
of optimizing the lattice constant is consistent with other 
approximations used in the calculations. 



B. DFT with vdW-DF 

In order to account for the van der Waals interaction 
and the effects of the long-range correlation in general, 
we applied the vdW-DF approach in a post-GGA proce- 
dure to DFT results for the structures (a), (c) and (d). 
Here the full power of the vdW-DF correlation approach 
becomes obvious, because due to its "seamless" character 
we can apply it at all graphene-Ir distances, i.e. at all 
coupling strengths, without worrying that it may spoil 
the GGA results which are good for strong bonding^ 

The results for structure (a) are shown by squares in 




FIG. 6: Standard DFT (LDA and PBE GGA) and vdW-DF 
energies for graphene monolayer on Ir(lll), structure (a) in 
Fig. [2 



> 

A 



Strong bond 



vdW minimum 
PBE 



vdW-DF 



LDA 



- Region of a lenthening 



> 

& 



° Strong bond 




vdW minimum 
PBE 



vdW-DF 



LDA 



Region of fl Q lengthening 



FIG. 7: Standard DFT (LDA and PBE GGA) and vdW- 
DF energies for Ir-graphene-Ir sandwiches, structures (c) and 
(d) in Fig. [4] The lines connect points for which the energy 
has minimum when allowing the structures to slightly expand 
laterally, as explained in the text. 



Fig. [6l where the PBE and the LDA results are the same 
as in Fig. 03 and the energy calculated using vdW-DF 
is shown by a thick line and squares. A clear van der 
Waals attractive well develops, deeper than the shallow 
well in the LDA calculation, and with the minimum at a 
larger graphene-substrate distance, around 3.7 A. Fig. [7] 
shows similar results for structures (c) and (d). The 
unconnected points arc for the lattice constant of free 
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graphene, as in Fig. and the points connected by lines 
for the optimized expanded lattice constant. The van 
der Waals potential well is similar to Fig. [5] but approxi- 
mately two times deeper (note the different scale on the 
energy axis), since the graphene interacts both with the 
iridium substrate and overlayer. The depth and shape of 
the chemisorption minimum at around 2.3 A is less af- 
fected, but the barrier between the two minima is much 
decreased compared with the DFT results. 



C. Discussion of the results 

Detailed information about the nature of C-C and C- 
Ir bonds can be inferred by examining the geometry of 
the graphene lattice around the minima of the interac- 
tion energy in Fig. [Jj At the physisorption minimum, 
Zh-g <~ 3.7 A, the graphene lattice is perfectly planar, 
and the graphene stays at midpoint between two neigh- 
boring iridium layers. The same is true for all struc- 
tures in Fig. [4j and in fact for other geometries such 
as monolayer graphene over the 1-2 and 2-3 regions in 
Fig. [3J Due to the smoothness of the potential with re- 
spect to the translation of graphene along the surface, 
graphene flakes physisorbed on Ir(lll) are quite mobile, 
both translationally and rotationally, which is an im- 
portant mechanism in aggregation and growth of large 
graphene islands^! 
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TABLE I: Bond length and angles at the strong bonding 
energy minima of structures (c) and (d). All lengths are in 
A. Here ao is the optimal lattice constant of the structure in 
the x — y plane, slightly larger than the free graphene value of 
3.478 A, zs-a is the distance between the substrate atom Irs 
and Ca, 2buck is the buckling of graphene, i.e. the difference of 
z coordinates of atoms Ca and Cb, zb-o the distance between 
Cb and the overlayer atom Iro, cLa-b the distance between 
two neighboring C atoms, and a the angle defined by the lines 
Irs-CA and Ca-Cb- 

The situation is rather different when strong bonding 
between iridium and graphene in Ir-graphene-Ir struc- 
tures occurs, at zi r _ g around 2.3 A in Fig. [7j First, we 
notice that the total energy depends strongly on the po- 
sition of the C atoms of graphene with respect to the Ir 
atoms below and above. Thus the two similar structures, 
(b) and (d) in Fig. [4] differ only in the position of the 
iridium atoms in the monoatomic overlayer, but the total 
energies in Fig. [5] differ by more than 1 cV! The absence 
of a stable strong bond in structure (b) shows that strong 
bonding can occur only when both C atoms are saturated 
by Ir atoms, one directly below and one above it. This 
immediately implies that the onset of strong bonding ef- 
fectively anchors the iridium cluster and the underlying 



graphene to the particular spot of the moire pattern of 
the graphene-substrate supercell. 

The formation of the strong "organomctallic" bond 
is accompanied by buckling of graphene and C-C bond 
lengthening. Table |T] shows the values of bond lengths 
and angles corresponding to the strong bonding energy 
minima in the two panels of Fig. [Jj These values are 
close to the ones of the tctrahedrally bonded C atoms in 
diamond, and indicate that the rehybridization from sp 2 
to sp 3 bonding has occurred, as noted by FeibelmanjS 



V. DISCUSSION 

These results show that the onset of the strong C- 
Ir binding in Ir-graphene-Ir sandwiches is accompanied 
by the disappearance of the aromatic character of the 
carbon rings. The carbon atoms rehybridizc to sp 3 con- 
figuration, and the two C atoms in the graphene unit 
cell move out of the plane in opposite directions. On 
the other hand, in one-sided binding on Ir(lll) (i.e. a 
clean graphene overlayer) the C-Ir bond is always weak, 
dominated by van der Waals interaction. 

In order to get more insight into the character of the 
bonding of aromatic rings on Ir(lll), we have made DFT 
calculations of benzene molecules CgHg lying flat on the 
iridium surface. We found large differences of binding 
energies and distances for different positions of the ben- 
zene molecules with respect to the substrate atoms. The 
most stable configuration is when the center of the ring is 
above a hollow site, and the six C atoms are above three 
Ir atoms, two C on each Ir. (This configuration cannot 
be directly compared to any part of the moire pattern 
of graphene on Ir(lll), Fig. G2 since it corresponds to a 
different orientation of the aromatic rings, i.e. rotated 
by 30°.) The C-Ir bond is around 2.4 A. The GGA ad- 
sorption energy is somewhat less than 1 cV, and does not 
change substantially when the vdW-DF nonlocal corre- 
lation is used. The C atoms remain planar due to sym- 
metry, but the C-C bonds become longer, 1.43 A and 
1.48 A, and the H atoms are slightly above the plane 
of the C atoms. This indicates a change of the nature 
of the bonding of the carbon ring and a departure from 
the pure sp 2 hybridization. The configurations where 
the center of the benzene ring is above an Ir atom, and 
H atoms point either towards the neighboring Ir atoms 
or towards bridge sites, are more weakly bound. The 
C-Ir bond length is around 3.4 A. In this case the non- 
local correlation is essential for the bonding. With pure 
GGA functional there is virtually no bonding of the ben- 
zene molecule, while with the vdW-DF the bonding en- 
ergy is around 0.6 eV. The C-C bonds keep the value of 
1.41 A as in the benzene molecule, and the whole benzene 
structure is planar. These values are also very similar to 
those of graphene on Ir(lll) obtained earlier. All these 
results indicate a weak, wan der Waals-dominatcd bond- 
ing. Thus the bonding of benzene on Ir(lll) shows even 
more variation of bonding parameters than various re- 
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gions of the moire of graphcne, indicating the richness of 
possible bonding of aromatic structures (molecules and 
graphene) on metal surfaces. 

Graphcne strongly binds on some other surfaces, 
as mentioned in the Introduction, apparently without 
strong rchybridization to sp 3 . In a recent combined ex- 
perimental and theoretical study of graphene bonding on 
Ru(0001) surface^ it was found that graphene is strongly 
corrugated with a minimum C-Ru distance of 2.1 A and 
a corrugation of 1.53 A in the regions of strong coupling. 
The DFT calculations were performed using the stan- 
dard PBE functional, which is expected to work well in 
the regions of strong coupling, and the lack of the van 
der Waals interaction which should be dominant in the 
weakly coupled 'blisters' is not crucial. The authors find 
that the height difference between neighboring C atoms 
in the graphene layer is below 0.03 A in the strong cou- 
pling region in DFT, and conclude that the adsorbed 
graphene layer remains sp 2 hybridized. 

Returning to the calculations of graphene on Ir(lll), 
we note that the use of the large supercell in Ref . HI has the 
advantage that the lattice constants of both the iridium 
substrate and the graphene overlayer can be kept close 
to their natural values. Thus the problems which we 
encounter with our compressed Ir surface (expanded in 
the z direction) are avoided. In particular, it seems that 
we get somewhat too large Ir-graphene distances com- 
pared to other calculations and the preliminary exper- 
imental estimates. Furthermore, when iridium clusters 
are added on top of graphene, in the large supercell ap- 
proach the carbon atoms are free to relax both vertically 
and laterally, which is essential for a good description of 
graphene buckling and the formation of the strong C-Ir 
bond. In subsection IIVBI we had to use a calculational 
tour de force to detect the preference of carbon atoms 
to lengthen somewhat the C-C bonds and thus approach 
more closely the diamond structure. Furthermore, in the 
supercell approach the substrate iridium atoms may also 
relax laterally, optimizing the saturation of the bonds to 
carbon atoms. 

These advantages come with the downside that in 
Ref. HI the LDA functional was used. This was a nec- 
essary choice since GGA in the usual formulation, i.e. 
without the van der Waals interaction being somehow 
accounted for, gives little or no binding of graphene (see 
Ref. HI, in particular the Erratum). However, LDA usu- 
ally gives a too small equilibrium distance, and overbinds 
in cases of strong chemical bonding. Thus in our calcu- 
lations of graphite in section [TH Fig. [2j LDA gave a too 
small intcrlayer distance. The LDA binding energy of 
graphite was also too small since graphite is a system 
with very little chemical component of the bond and the 
LDA ovcrbinding could not compensate fully the absence 
of the vdW component. 

The compressed Ir(lll) surface in our approach and 
the use of LDA in Ref. IH preclude a detailed quantitative 
comparison between the results of the two calculations, 
and of each of them with experiment. However, the semi- 



quantitative agreement is good. Both approaches predict 
a rather weak bonding of a graphene monolayer with the 
Ir(lll) substrate, and the formation of a much stronger 
organometallic bond when iridium clusters are added 
on top, accompanied with the buckling of the graphene 
structure and shortening of Ir-C distances. For clean 
graphene, the Ir-C separation at the 1-3 regions of the 
moire pattern is around 3.48 A in Ref. |9| and around 
3.7 A in our work, while the experimental value has been 
estimated to around 3.38 A. When the clusters trigger 
strong bonding and graphene buckling the Ir-C distance 
decreases to around 2.1 A in Ref. and to around 2.2 A 
in our work, while the Ir-C-C angles are around 105°. 

The bonding of graphene on some other (111) surfaces 
of fee metals assuming commensurate configurations has 
also been investigated. In the papers by Giovannetti et 
al.— and Khomyakov et al.,— the LDA functional was 
used. The unit cells were either 2 graphcne C atoms and 
one metal atom in each layer (e.g. Ni, Co, Cu), as in 
our calculation, or 8 C atoms and 3 metal atoms with 
the graphene unit cell rotated by 30° when the difference 
of the lattice constants was larger (e.g. Pd, Au, Pt). 
It was found that graphcne interacts strongly with Ni, 
Co, and Pd, with the equilibrium metal-graphene dis- 
tance between 2.05 and 2.30 A, and weakly with Cu, 
Au and Pt, with equilibrium distance between 3.26 and 
3.31 A. These findings are in agreement with experiment, 
where available. The mismatch of the lattice constant of 
graphene is 4% for Cu, 1.2% for Ni and 2% for Co (metal 
unit cell being larger in all three cases). This is clearly 
smaller than in our calculation, where the difference of 
Ir(lll) and graphene lattice constants is around 10%. 

Vanin et ali22 consider the same surfaces, but in a 
quite different approach. They use the vdW-DF correla- 
tion functional^ evaluated using the method proposed in 
Ref. and self-consistently implemented into the real- 
space projector augmented wave gpaw codeJ^ They do 
not adjust the metal substrate to match the lattice con- 
stant of graphene, but instead keep it at their experi- 
mental lattice parameters and adjust the graphene sheet. 
They claim that the vdW-DF results do not change sig- 
nificantly if they fix the graphene lattice parameter to its 
optimized value and adjust the metals correspondingly. 
Surprisingly, they obtain weak binding for all metals con- 
sidered, with metal-graphene distances in the range 3.40- 
3.72 A. This is in clear disagreement for Co and Ni, where 
strong binding has been experimentally confirmed. 

In contrast to this, our calculations of Ir(lll)-graphcnc 
structures give an overall agreement with other calcula- 
tions and with experiment. The weak vdW minimum 
around 3.7 A which exist both for clean graphene over- 
layer (Fig. [6]) and for graphene with iridium adclusters 
(Fig. [7|) is obtained correctly only with vdW-DF. The 
overall shape of the potential minimum of the strong 
bond around 2.3 A for iridium adclusters is roughly sim- 
ilar for calculations using LDA and PBE with vdW-DF, 
although LDA clearly overbinds. Even plain PBE calcu- 
lations show a comparable local minimum, just weaker. 
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The disagreement found in Ref. l22j is therefore even more 
surprising. We have not tried our method on metals other 
than iridium, where the strong chemisorption minimum 
exists only if adclusters are present. There is a possibil- 
ity that vdW-DF does not work so well for other metals. 
However, in our opinion the source of disagreement may 
also be the other approximations used in Ref. [22T 

The graphenc is particularly stable due to the aro- 
matic character of the carbon rings. Perturbing the 
structure (for example by forcibly changing the natu- 
ral bond length) may significantly change the reactivity. 
Thus simply adapting the graphene lattice constant to 
the substrate may have unwanted consequences, weaken- 
ing the stability of the aromatic bonds, as well as chang- 
ing the doping of the graphene layer in contact with the 
metal surface. The opposite procedure, which we used in 
this paper, i.e. adapting the substrate lattice constant, 
seems preferable to us but may also have some weak- 
nesses. First, the change of the electronic structure of 
the substrate may be large enough to alter the reactivity 
compared to the natural metal. Also, the lattice constant 
of the free graphene may not be optimal for rehybridized 
graphenc forming strong sp 3 bonds. We had to expand 
the graphene lattice slightly in order to obtain a suffi- 
ciently stable strong bonding. In the process we had to 
carefully account for the change in energy of the iridium 
substrate, which was, of course, also expanded (i.e. less 
compressed compared to the natural structure) . All this 
indicates that the graphene lattice constant should be 
left at its natural value in the weak bonding cases, but 
should be allowed to relax and lengthen when the strong 
bonding regime accompanied with graphene buckling and 
rehybridization to diamondlike bond occurs. This cannot 
be achieved in the simplified commensurable geometries, 
and a full large supercell calculation with state-of-the-art 



nonlocal correlation functional seems to be the only ap- 
proach which can give the answer about the structure of 
graphenc adsorbed on various metals in the general case. 



VI. CONCLUSIONS 

We find that a graphene monolayer on Ir(lll) is weakly 
bound, and keeps the aromatic character of the car- 
bon rings. In Ir-graphene-Ir structures C atoms show 
a tendency towards rehybridization and formation of sp 3 
bonds, which in favorable cases (an Ir atom directly be- 
low or above each C atom) are more stable that the ph- 
ysisorbed structure. In all cases, the use of the vdW-DF 
which includes a full description of the nonlocal corre- 
lation is essential. However, our approach in which the 
substrate lattice constant is adjusted to match graphenc 
does not give full quantitative accuracy. In order to ob- 
tain that kind of agreement, large calculations on realis- 
tic supercells using DFT functionals with nonlocal cor- 
relation are necessary. This conclusion is also true for 
other graphene-on metal systems, in which the nature of 
the graphene-metal bond may be quite different than on 
Ir(lll). 
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